{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "9d9f3884-38ee-412b-80cc-71bd6c86276f",
   "metadata": {},
   "source": [
    "# The `Reiter` notebook  <a id=\"Reiter\"></a>[<font size=1>(back to `Main.ipynb`)</font>](./Main.ipynb)\n",
    "\n",
    "This notebook gathers the functions for the comparison with the Reiter method.  \n",
    "\n",
    "The function \n",
    "\n",
    "> `Reiter(solution::AiyagariSolution, economy::Economy)`\n",
    "\n",
    "returns a tuple spcific to the comparison with Reiter.\n",
    "\n",
    "The function \n",
    "\n",
    "> `write_Reiter(filename::String, solution::AiyagariSolution, economy::Economy)`\n",
    "\n",
    "writes the output in a `mat` file `filename` for the comparison with Reiter. It calls the `Reiter` function to compute the specific elements.\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "id": "54cd646f",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Reiter (generic function with 1 method)"
      ]
     },
     "execution_count": 99,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function Reiter(solution::AiyagariSolution, economy::Economy)\n",
    "    @unpack β,α,δ,χ,φ,σ,Tt,u,u′,na,a_min,aGrid,ny,ys,Πy,l_supply = economy\n",
    "    @unpack ga,gc,R,w,A,Y,K,C,L,transitMat = solution\n",
    "    I = typeof(na)\n",
    "    T = typeof(β)\n",
    "            \n",
    "    Va = ga[:]'\n",
    "    Vind = zeros(Int,na*ny)\n",
    "    Wind = zeros(Float64,na*ny)\n",
    "\n",
    "    for (i,a) in enumerate(Va)\n",
    "        set = findall(x->x<=a,aGrid)\n",
    "        if isempty(set)\n",
    "            Vind[i] =1 \n",
    "        else\n",
    "            Vind[i] = min(maximum(findall(x->x<=a,aGrid)),length(aGrid)-1)\n",
    "        end\n",
    "        Wind[i] = 1 - (a - aGrid[Vind[i]])/(aGrid[Vind[i]+1]-aGrid[Vind[i]])        \n",
    "    end\n",
    "\n",
    "    gc_Reiter =  gc[:]\n",
    "    lv = kron(l_supply.(w,ys),ones(na,1))\n",
    "    \n",
    "    x = gc[:] - (1/(χ*(1/φ+ 1)))*(lv).^(1+1/φ)\n",
    "    resE = zeros(Float64,na*ny)\n",
    "    for i=1:na\n",
    "        for j=1:ny\n",
    "            ind = i+(j-1)*na\n",
    "            EUc = zero(Float64)\n",
    "            for jp=1:ny\n",
    "                EUc += Πy[j,jp]*(\n",
    "                        Wind[ind]*x[Vind[ind]+(jp-1)*na] +\n",
    "                        (1-Wind[ind])*x[Vind[ind]+1+(jp-1)*na])^-σ\n",
    "            end\n",
    "            resE[ind] = x[ind]^-σ- β*R*EUc\n",
    "        end\n",
    "    end\n",
    " \n",
    "    # Computing the consumption share of the bottom 10%\n",
    "    M_temp = [repeat(aGrid,1,ny)[:]  solution.stationaryDist[:]]\n",
    "    indsort = sortperm(view.(Ref(M_temp), 1:size(M_temp,1), :))\n",
    "    cumDist = cumsum((solution.stationaryDist[:])[indsort]) # cumulative stationary distribution \n",
    "                                                            # sorted by beg-of-period wealth 1st and \n",
    "                                                            # share in the population then.\n",
    "    ind_bot10 = indsort[findall(x -> x<=0.10, cumDist)] # indices of the bottom 10%\n",
    "    c_bot10 = sum(solution.stationaryDist[ind_bot10].*gc[ind_bot10])/sum(solution.stationaryDist[ind_bot10])    \n",
    "                    # per-capita consumption of the bottom 10%\n",
    "    \n",
    "    return gc_Reiter,Vind,Wind,resE,Va,ind_bot10,c_bot10\n",
    "end;"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 106,
   "id": "66f6a408",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "write_Reiter (generic function with 1 method)"
      ]
     },
     "execution_count": 106,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "function write_Reiter(filename::String, solution::AiyagariSolution, economy::Economy)\n",
    "    gc_Reiter,Vind,Wind,resE,Va,ind_bot10,c_bot10 = Reiter(solution,economy) \n",
    "    \n",
    "    economy_Dynare = Economy_Dynare(economy)\n",
    "\n",
    "    file = matopen(filename, \"w\")\n",
    "    write(file, \"eco\", economy_Dynare)\n",
    "    write(file, \"solution\", solution)\n",
    "    write(file, \"Vind\", Vind)\n",
    "    write(file, \"Wind\", Wind)\n",
    "    write(file, \"resE\", resE)\n",
    "    write(file, \"Cpol\",gc_Reiter)\n",
    "    write(file, \"b\",ind_bot10)\n",
    "    write(file, \"Cbs\",c_bot10)\n",
    "    close(file)\n",
    "    return nothing    \n",
    "end;"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Julia 1.10.5",
   "language": "julia",
   "name": "julia-1.10"
  },
  "language_info": {
   "file_extension": ".jl",
   "mimetype": "application/julia",
   "name": "julia",
   "version": "1.10.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
